Faddeev calculations of break-up reactions with realistic 

experimental constraints 

J. Kuros-Zolnierczuk^, P. Thorngren Engblom^, H-0. Meyer\ T. J. Whitaker^ 
H. Witak^, J. Golak^, H. Kamada^ A. Nogga^ R. Skibinski^ 

^Department of Physics, Indiana University, Bloomington, IN ^l^OS, USA 
'^Department of Radiation Sciences, 
Uppsala University, S - 75121 Uppsala, Sweden 
^M. Smoluchowski Institute of Physics, 
Jagiellonian University, Reymonta 4, 30-059 Krakow, Poland 
'^Department of Physics, Faculty of Engineering, 
Kyushu Institute of Technology, Kitakyushu 804-8550, Japan and 
^ Institute for Nuclear Theory, University of Washington, 
Box 351550, Seattle, WA 98195, USA 

Abstract 

We present a method to integrate predictions from a theoretical model of a reaction with three 
bodies in the final state over the region of phase space covered by a given experiment. The 
method takes into account the true experimental acceptance, as well as variations of detector 
efficiency, and eliminates the need for a Monte-Carlo simulation of the detector setup. The method 
is applicable to kinematically complete experiments. Examples for the use of this method include 
several polarization observables in dp break-up at 270 MeV. The calculations are carried out in the 
Faddeev framework with the CD Bonn nucleon-nucleon interaction, with or without the inclusion 
of an additional three-nucleon force. 
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I. INTRODUCTION 



For nuclear reactions with two particles of given masses in the exit channel, two parame- 
ters are required to specify the final state. Usually, these are taken to be the polar angle and 
the azimuth of one of the particles. The direction of the other particle and both energies 
are then given by energy and momentum conservation. 

When three particles are present in the final state, one more momentum vector needs to 
be specified, thus three more parameters are required. This means that all observables are a 
function of five parameters. Often these are taken to be the polar angles and the azimuths 
of two of the three particles and one parameter that describes how the kinetic energy is 
shared by the three particles. This parameter could be, e.g., the relative energy of particles 
1 and 2. 

In order to present and discuss experimental results, one must select one or two inde- 
pendent variables. Traditionally, two outgoing particles are measured in coincidence by two 
(small) detectors. The detector positions determine the polar and azimuthal angles of both 
particles. In a plot of the energy of the first versus the energy of the second particle, all 
events lie on a locus. The position along that locus then serves as the independent variable. 
The choice of the detector angles is arbitrary, but their solid angles must be reasonably small 
to keep the locus defined. To explore the full kinematics of the reaction the experiment must 
be repeated with different detector positions (see e.g. 

Modern nuclear detection techniques make it possible to build experiments that simul- 
taneously cover a sizable fraction of the entire phase space of the reaction. Measuring the 
momenta of two outgoing particles yields enough information to define the complete kine- 
matics of an event, with one redundant parameter to spare that may be used for event 
identification. One thus has the freedom of choosing any independent variable by sorting 
the measured events accordingly. 

When partitioning the phase space with respect to a given single variable, all other 
kinematic parameters are ignored. When one wants to compare the experiment to a model, 
the corresponding calculation must integrate over these ignored parameters. To carry out 
this integral, one must take into account the boundaries of the acceptance of the experiment, 
as well as variations of the detection efficiency within the acceptance. It is usually quite 
difficult to describe the detailed performance of the detection system, and to implement this 
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information in the calculation of an observable by a theoretical method. In this paper we 
present a straight-forward method to accomplish this task. 

In Sec.|n]we describe the proposed method. In Sec. Illll we outline the underlying theory 
of three-nucleon break-up reactions and discuss how the calculations can be accelerated for 
the use in the present context. In Sec. IIVI the new method is applied to the analysis of 
several polarization observables in the dp break-up reaction. These calculations are based 
on the charge-dependent (CD) Bonn nucleon-nucleon (NN) potential . We also study 
the effect of including the modified Tucson-Melbourne three-nucleon force (TM' 3NF) |5[ . 
Finally, we summarize our results in Sec. |^ 



II. THE SAMPLING METHOD 



Let us denote by x = {ai, •••} the set of parameters that is needed to completely 

describe the kinematics of a given nuclear reaction. To specify the phase-space coordinates 
of a three-body-final state one requires 5 parameters. 

The differential cross section for the reaction with unpolarized collision partners is given 
by a"o(x). A typical polarization observable 0{x) has the effect of modifying the unpolarized 
cross section such that (7(x) = cro(a;)(l + PO{x)), where P is the polarization of the beam or 
the target or their product, and 0{x) is a beam analyzing power, a target analyzing power 
or a spin correlation coefficient. In order to measure 0{x) one carries out two measure- 
ments with opposite sign of the polarization P. The yields N± accumulated during the two 
measurements in a phase space region [x, x + Ax) are then given by 

N±{x) = L±e{x)ao{x){l ± P±0{x)) (1) 

Here, P+ and P_ are the magnitudes of the two polarizations with opposite sign, and L_|_ and 
L__ are the time-integrated luminosities for the two measurements. The detector efficiency 
e{x) measures the probability with which an event of interest gets registered by the detector 
system. For events outside the acceptance of the detector, e(x) =0. 

Let us assume, for the moment, that the integrated luminosities and polarizations for the 
two measurements are the same, or L+ = L_ = L/2, and P^ = P^ = P. We then find for 
the total number of collected events 

N{x) = N+ + N_ = Le{x)ao{x) (2) 
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and for the polarization observable in terms of the measured yields 



0(x) = (l/P)(iV+ - iV_)/(iV+ + N.) (3) 

The equation for O holds for any point x in phase space, but it is obviously impractical 
to evaluate or discuss observables in a five- dimensional parameter space. Rather, one may 
select a single independent parameter am, ignoring all others. When the experiment ignores a 
given parameter, the full range of that parameter within the detector acceptance is included. 
Thus, for each of the Om bins, one actually evaluates 0(7) in a region 7 of phase space, 
where 7 denotes the range for each of the ignored parameters. Assume now that we wish 
to compare an experiment 0(7) to a theoretical model. The model calculation provides 
us with a value 0^^{x) at any point x in phase space. In order to obtain the theoretical 
equivalent 0*^(7) of the experiment we have to average 0*'^{x), over the region 7, weighted 
by the unpolarized cross section times the experimental efficiency, 

f ao(x)e(x)0'^'^(x)dx 

0'\l) = ' ' ■ (4) 

J^Mx)(^{x)dx 

We determine e(x)(Jo(x) by making use of Eq. (j2)), and we replace the integrals in Eq. (jl} 
by sums over all elements Xi of size Ax, that make up the region 7, 

Here, N{xi) is the number of events collected in element Xi, irrespective of polarization. 
Since we are free to choose the size Ax of the element, we decrease it until all N{xi) are 
either or 1. The number of occupied x elements then equals the total number of events 
A^(7) collected in region 7 during the experiment, and the list of the Xj's for occupied bins 
is then identical to the list of phase space coordinates Xk {k = l..N{'y)) for all collected 
events. For a kinematically complete experiment, in which the phase space coordinates x^ 
are known for each event, the correctly averaged value for a calculation is thus obtained as 
the mean of the corresponding theoretical values O*'' for all events k 

0'\^) =< O'^ >= (6) 

This simple recipe constitutes our proposed method: in order to obtain the average theo- 
retical value, correctly weighted by the product of unpolarized cross section and detector 
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efficiency, one determines the theoretical O for each collected event, sums these values and 
divides by the number of events. 

It is easy to see, that the standard deviation of the theoretical value that arises from the 
randomness of the experimental phase space points that are used to sample the region 7 is 
given by 



The method of Eq. 0, in effect, involves an average over theoretical values, weighted 
by the density in phase space of the actual events. This density is proportional to the 
unpolarized cross section and detector efficiency, but it also depends on the polarization 
asymmetry. The latter cancels only if the magnitudes of the two polarizations with opposite 
sign are the same, and the corresponding, time-integrated luminosities L_|_ and L_ are the 
same, or AP = (P+ - P_)/2 and AL = (L+ - L_)/2 both vanish. If AP ^ 0, then 
the averages < O*^ >^ and < O*'* >~, taken with just the data points with positive or 
negative polarization, respectively, are different. It is easy to see that in this case, the 
desired theoretical average < O*^ >o that is free of polarization effects can be obtained 
from the average < O*'* > obtained with all events according to Eq. ©, by subtracting a 
correction term, 



In the experiment [6| to which we later apply our method, an effort was made to keep AP 
and AL small, and the correction term of Eq. |H1 turned out to be insignificant. 

III. THEORY USED TO CALCULATE OBSERVABLES 
A. Basic definitions 

As an example of an application of the method described in this paper, let us consider 
the proton-deuteron break-up reaction at an energy below the pion production threshold. 
Here, all three nucleons are moving freely in the outgoing channel. The kinematics of a three 
body final state is determined by nine variables, but energy and momentum conservation 
reduces the number of independent variables to five. 
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Assume that a.d + p-^p + p + n is kinematically complete, as is the case when the 
energies and directions of the two final-state protons are detected. To describe the final 
state kinematics we define the Jacobi momenta p and q: 



where bi and 62 are momenta of the two protons in the center-of-mass system (see Fig. QJ. 



FIG. 1: The momenta of three particles in center-of-mass system. Particles numbered 1 and 2 are 
the two protons with momentum bi and 62- 

In the remainder of this paper we are concerned only with 'axial' observables that are 
symmetric with respect to a rotation around the beam axis. For this special case, the number 
of independent variables reduces to four, since only the difference between the azimuths is 
relevant. We choose these variables to be a; = {p,9p,9q, Acj) = (pp — 0^}, in the center of 
mass. 

B. The break-up amplitude 

The theoretical predictions presented in this work are based on solutions of 3N Faddeev 
equations using realistic NN interaction and a 27r-exchange three-nucleon force model. In 
the following, we give a short overview of the underlying formalism. For more details we 
refer to 0, 0| and references therein. 

The transition operator for the break-up process Uo can be expressed in terms of a T 
amplitude as: 




q = -{bi + h) 



(9) 





(10) 
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The operator T fulfills Faddeev-like equation 

T = tV + {l+ tGo)v}^\l + V)+ tVGoT + (1 + tGo)v}^\l + V)GqT (11) 

where the two-body t-operator is denoted by t, the free 3N propagator by Go and V is the 
sum of a cyclical and an anti-cyclical permutation of three particles. The three-nucleon force 
V4 can always be decomposed into a sum of three parts: 

V, = V^'^ + Vt^ + Vt\ (12) 

where V^^*^ singles out nucleon i and is symmetrical under the exchange of the other two 
ones. As seen in Eq. flll|) only one of the three parts occurs explicitly, the others enters via 
the permutations contained in V. 

As a 2N-force we use the CD Bonn NN potential This interaction is one of the 

modern, phenomenological models which describe the 2N data set with a per data point 
that is very close to one. As a 3NF we take the new version of the Tucson-Melbourne force - 
TM' 3NF 0,0]. The strong cutoff parameter A has been adjusted such that the 3NF together 
with CD Bonn potential reproduces the measured triton binding energy (A = 4.593m7r jsj). 

We solve Eq. (jllj) in a partial-wave projected momentum-space basis 0|. The calcula- 
tions shown here are carried out for a deuteron bombarding energy of Efi =270 MeV (the 
equivalent laboratory energy of incoming proton is Ep = 135 MeV). In order to achieve 
converged solutions of the Faddeev equations, it is necessary to include partial waves up to 
the 2N-subsystem total angular momentum jmax = 5. This corresponds to up to 142 partial 
wave states in the 3N system. The 3NF includes all total angular momenta of the 3N system 
up to J = 13/2 while the longer-ranged 2N interactions require states up to J = 25/2 for 
convergence. The Coulomb force between the protons is neglected in these calculations. 
Coulomb effects are expected to be small at the considered energy of the incoming proton 
and to show up only in final-state configurations with low relative proton momenta. 

The Faddeev calculation results in the break-up amplitude Uo- From this amplitude 
one obtains, in a standard manner, the cross section or any polarization observable [ij. 
Definitions of polarization observables with a spin-1 and a spin-1/2 particle in the initial 
state can be found in j^. 

In the following, we will discuss the longitudinal proton analyzing power Az which is 
non-zero for non-coplanar p and q, but forbidden by parity conservation when equals 



or 71. In addition we have calculated the deuteron tensor analyzing power Azz and some spin 
correlation coefficients Ca,b, that can be measured with both initial-state particles polarized. 
Here, a refers to the vector or tensor polarization of the deuteron, and b to the polarization 
direction of the nucleon. 

C. Computational details 

In Sec. |n] we have presented a method to calculate the theoretical estimate for any 
observable, correctly weighted by the experimental acceptance and detector efficiency. The 
method requires that a theoretical value is obtained for the phase space coordinates of 
every measured event. However, deducing such a value directly from the Faddeev amplitude 
(single-shot mode) for typically a few million events is not practical. In order to overcome 
this limitation, we precalculate and store the desired observable (0*''(x)) at discrete points, 
covering the entire phase space. The value of the observable at any phase space point is 
then retrieved by multidimensional linear interpolation as described in the Appendix. 

The observable is precalculated at uniformly-spaced phase-space points, called a 'grid', 
chosen as follows. The azimuth variable A0 is taken in steps of 10° between and 360°. 
The polar angles 6p and 6q are taken in 5° steps, from 5° to 90° and 180°, respectively. For 
6p only half the angular range is needed because the two observed particles in the final state 
are identical. The range of allowed p values is divided into 30 steps. Thus, the grid points 
form a matrix with 719,280 elements. 

The required mesh size has been determined such that the shape and magnitude of 
< O*^ > is sufficiently insensitive to a change of the number of grid points. In order to test 
how the interpolated values depend on the mesh size of the grid, we performed calculations 
using grids with different numbers of points. The average of Eq. (jH)) is always taken over the 
same list of about 3 * 10^ phase space coordinates, culled from an actual experiment 0|. In 
Fig. 121 we compare the analyzing power < A*^^ > calculated with grids of fewer points (coarser 
mesh size) to the corresponding < A*'* >^*^ obtained with the fine grid described above. The 
ratio of the two values is shown as a point for each of the 36 bins of the independent variable 
A(f). A vertical spread of points is thus a measure for the interpolation uncertainty. 

To compare the (time-consuming) single-shot mode with the (fast) grid interpolation 
method, we have calculated the analyzing power < A'^J^ > as a function of A(f) using both 
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FIG. 2: The ratio between the analyzing power < ^* > calculated with grids with smaller numbers 
of grid points and the analyzing power < A''J^ >^*^ obtained with the final, fine grid. The ratio of 
the two values is shown as a point for each of the 36 bins of the independent variable Ac/). 

methods. The values of A^J^ for all events within a given A0 bin are averaged according to 
Eq. (jUj) to give < Az > for that bin. The error bars are calculated according to Eq. ((7j) 
and turn out to be smaller than 0.003. The results from the single-shot mode (circles) and 
the grid interpolation method (crosses) are shown in Fig. The difference between the 
two results is shown in Fig. ^jp. Again, this difference is a measure of the interpolation 
uncertainty, or, in other words, of the 'curvature' of the interpolated function (hence the 
systematic dependence on A^). 



IV. EXAMPLES 



Recently, an experimental study of the dp break-up reaction with polarized 270 MeV 
deuterons on a polarized proton target has been completed using the "Cooler" storage ring 
at the Indiana University Cyclotron Facility 0]. In this section, we illustrate the use of our 
new method to obtain theoretical estimates of polarization observables for this experimental 
setup. 

We calculated the differential cross section and several polarization observables 0^^{x) 
using the CD Bonn NN potential with and without the TM' 3NF for each grid point in 
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FIG. 3: a) comparison of the single-shot calculation of A^^ (open circles) with the grid interpolation 
method A™* (crosses); b) the difference between the two calculations 

the available phase space. After that, using the multidimensional interpolation, cf. the 
appendix, we obtained theoretical predictions for each experimental phase space point (about 
3 * 10^ experimental points). Using the averaging method of Eq. (jHI), the boundaries of the 
experimental acceptance for the energies and azimuths of the two outgoing protons, which are 
approximately E{^^,El^^ > 50 MeV, 10° < e^f,9^f < 45°, respectively, were automatically 
taken into account. The errors from the sampling statistics are obtained from Eq. ((7|). 

The results for the vector analyzing power A^, the tensor analyzing power A^^ and the 
correlation coefficients Czz,zi Cy^ — Cxy and C^x + Cyy are shown in Figs (HJ - © as a function 
of A0. In these figures the solid line shows the Faddeev calculations based on the CD Bonn 
NN potential, while the dashed line is obtained by including the TM' 3NF. The errors (less 
than 0.003) are too small to be visible in figures. 

The difference between the solid and the dashed lines illustrates the effect of the TM' 3NF. 
As one can see, the 3NF effects are rather small in the present case where the observables are 
integrated over a large portion of the phase space. However, in the case of some correlation 
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FIG. 4: Average value of vector analyzing power as a function of Acp. The solid and dashed lines 
are the CD Bonn and CD Bonn + TM' 3NF predictions, respectively. 



0.2 r 



— CD Bonn 

-- CD Bonn + TM' 



0.1- 




90 



180 270 

A4> [deg] 



360 



FIG. 5: Average values of tensor analyzing power as a function of A(f). The solid and dashed lines 
are the CD Bonn and CD Bonn + TM' 3NF predictions, respectively. 

coefficients (Fig. ^ the magnitude of 3NF effects is sufficiently large to be distinguishable 
by a future experiment. 

It is conceivable that the smallness of the 3NF effects in some of the cases that are shown 
here is due to the fact the observable is averaged over a large portion of the phase space. 
In order to identify regions in phase space with strong 3NF effects, we have searched for 
the largest differences between observables predicted with the CD Bonn interaction either 
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FIG. 6: Average values of spin correlation coefficients as a function of A(j). The solid and dashed 
lines are the CD Bonn and CD Bonn + TM' 3NF predictions, respectively. 

by itself, or with the TM' 3NF included. The influence of the 3NF can be quantifled by a 
difference AO*'' at a given phase space point when the 3NF is switched on 

^Qth ^ \o"'{2N) - 0"'{3N)\ 

Here 0*^{2N), 0^^{3N) are predictions for the particular observable resulting when using 
the CD Bonn potential alone and including in addition the TM' 3NF, respectively. 

First, for each observable O^'^ we searched over the entire grid of the phase space points 
for the maximum value of the difference AO*'*. The results are presented in Table I. Second, 
since we are looking for the effects that should be experimentally accessible, we searched for 
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the maximum of AO only through the phase space points for which the differential cross 
section = a(2N) > 0.001- 



dfxiqdq — ^ \"^' I ^ " MeVsr^ " '^^^ resulting 'max\lS.O \ are shown in third column 
in Table I. In addition, in the brackets, labeled (*) and (**) of Table I, we list the number 
of grid points for which 0.05 < AO*'' < max[AO*^] and 0.05 < AO*^ < max[AO*^] with 



a(2iV) > 0.001 



MeVsr'^ ' 



respectively. As we can see in Table I the tensor analyzing power 



TABLE I: Maximum value of difference AO*'' = |0*''(2iV)-0*''(3Af)| for each observable O^^ . The 
second column, labeled maa;[AO*''], gives results without any restrictions. The third one, labeled 
ma3;[A0*''] + cr, results after applying cross section limitation ( cr(2iV) > 0.001 ). In the 

brackets, in the column labeled (*), there is a number of grid points for which 0.05 < AO*'* < 
max[AO**'']. In the column labeled (**) there are number of grid points for which 0.05 < AO*'' < 
maa;[AO*''] and a(2iV) > O.OOl^g^^. 
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or the correlation coefficient Cxx + Cyy are especially interesting, showing large 3NF effects, 
on the other hand the vector analyzing power A^ shows the least evidence of the action of 
the TM' 3NF. 

The next step in our study was to find the regions of the phase space with the largest value 
of AO*''. Let us investigate here, only for example, the dependence of y4^^(A0) on the addi- 
tional variables p, 9p and 9q. For this purpose, the four-dimensional phase space distributions 
of < A*ll > and AA*'^ are projected onto three planes: A(f) — p, A0 — 6p and A0 — 6^. In the 
case of AA*^ for each projection point the maximum value over the range of the remaining 
variables is determined, e.g. the A(j)—p projection shows max[Ay4*'j(p, 9p, 9q, A0)]A<^,p. Using 
graphs from the second column we can see which region of phase space exhibits a large 3NF 
effect. The figure shows, for example, that an interesting region would be 100° < Og < 180° 
with the full range of p and 9p included. Fig. |H1 shows the average value of < A*^ > after 
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A^^ih(2N) max[IA^^th(2N)-A^^ih(3N)l] 





FIG. 7: Two-dimensional projections of theoretical prediction for tensor analyzing power AJl- 
Columns show the average value of < A^Jl{2N) > and the maximum of difference \Af^{2N) - 
A^^{3N)\. The three rows show projections onto the A(p — p, Ac/) — Op and Acf) — 9q planes, 
respectively. 

integration over all values of p, Op and integration over part of 100° < 9g < 180°. For this 
example the 3NF effects are quite large. 



Studies such as this, together with earlier work published in 



may be used to identify 



regions of enhanced sensitivity the the 3NF in order to plan future experiments. 



V. SUMMARY 



We have presented a simple, straight-forward method to integrate over regions of phase 
space the prediction from a theoretical model. The method reflects the true experimental 
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FIG. 8: The average value of tensor analyzing power Azz as a function of Acf), in a region of 
phase space that is sensitive to the inclusion of a three-nucleon force (0.05 < p < 1.45,5° < Op < 
90°, 100° <eg< 180°). The sohd and dashed hues are the CD Bonn and CD Bonn + TM' 3NF 
predictions, respectively. 

acceptance, as well as variations of the detector efficiency. It is applicable to kinematically 
complete experiments with three bodies in the final state, and eliminates the need for a 
Monte-Carlo simulation of the detector setup. 

We have applied this method to several observables of dp break-up at 270 MeV, and 
studied the effect of including a three-nucleon force in the model. 

Acknowledgments 

The work has been supported by the U.S. National Science Foundation under Grant NO. 
PHY-0100348, the Swedish Research Council Dnr 629-2001-3868, DOE grants DE-FC02- 
01ER41187, DE-FG03-00ER41132 and by the Polish Committee for Scientific Research un- 
der Grant No. 2P03B00825. The authors would like thank to Professor W. Gloeckle for 
fruitful discussions and remarks. One of us (R.S.) would like to thank the Foundation for 
Polish Science for a financial support. The Faddeev calculations have been performed on 
the SVl and the CRAY T3E of the NIC in Juelich, Germany. 




15 



APPENDIX A: MULTIDIMENSIONAL LINEAR INTERPOLATION 



The sampling method described in section |nj requires that the theoretical value of a given 
observable is calculated for each of a large number of events, collected in an experiment. In 
order to achieve the necessary processing speed, precalculated theoretical values are stored 
in a four- dimensional matrix (the 'grid') as explained in Sec. IIII CI Multidimensional linear 
interpolation is then used to retrieve the corresponding value at any phase space point. 

We denote the set of four continuous independent variables as x = {p, 9p, 9q, A0} = {xi}; 
{i = 1..4). The corresponding discrete independent variables that constitute the grid are 
Xi,ki, ki = l..n, where n is the size of the grid in the ith dimension. It is understood that 
Xi,ki ^ Xi < Xi^ki+i for any value Xi of the set x at which a functional value is sought, i.e. the 
elements of the discrete grid ajj^fcr'^ariables are strictly increasing and the Xi coordinate at 
which interpolation is performed lies always inside an interval of grid points. Applying linear 
interpolation in four dimensions, means that the grid cell that contains the argument 
X of interest has 2^ = 16 corners. We are using the FORTRAN routine FINT 3] from the 
CERN program library. The interpolation is carried out as follows. 

The interpolation weights, Wi, are defined as 

Wi = "^^"^'-'^ (Al) 

Xi,ki+1 •^i,ki 

For brevity we also define Ui = 1 — Wi. The functional values at the corners of the grid cell 
that contains x are denoted by Um'-, fn = 1..16, and defined as follows, 

l/l = 0{Xi^ki, X2,ki, X3^ki, X4^ki) 
1/2 = 0{Xi^ki+l, X2,ki, X^^ki, Xi^ki) 
y3 = 0{xi^k,+l, X2,k,+1, Xs^ka X^m) 

1/4 = 0(xi,fc,,X2,fc,+i,a;3,fc,,a;4,fcJ 

1/5 = 0{Xi^ki, X2^ki, X3^ki+1, Xi^ki) 

ye = 0(xi,fc,+i,a;2,fc,,a;3,fc,+i,a;4,fcJ 

1/7 = 0(xi,fc,+i,X2,fci+l,X3,fc,+i,X4,fcJ 
Us = 0{Xi^ki, X2,ki+1, X3^ki+1, X4^ki) 



(A2) 



J/16 = 0{xi^k,, X2,k,+1, X3^k,+1, X4,,ki+l) 
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The observable at any point x = (xi, X2, X3, X4, X5) is then given by 

0{Xi,X2,Xs,X4) = 

U1U2U3U4 ■ yi + W1U2U3U4 ■ 1/2 + W1W2USU4 ■ ys + U1W2U3U4 ■ 1/4+ 
U1U2W3U4 ■ 2/5 + W1U2W3U4 ■ ye + W1W2W3U4 ■ yj + U1W2W3U4 ■ 1/8+ (A3) 
U1U2U3W4 ■ yg + W1U2U3W4 ■ yio + W1W2U3W4 ■ yii + U1W2U3W4 ■ yi2+ 
U1U2W3W4 ■ 2/13 + W1U2W3W4 ■ yi4 + W1W2W3W4 ■ 2/15 + U1W2W3W4 ■ yiQ 
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FIG. 9: The difference between the interpolated and the exact test observable. 

In order to investigate the accuracy of the interpolation routine, we construct a 'test' grid 
that is analogous to the one described in Sec. IIII C\ but with function values that are given 
by the analytic expression 

0{X) = A^_testiP, Op, Oq, 0p) = ^^^^ 

4:p/pmax{l -p/Pmax) " (sin(6'g) sin(26'p) sin(0p) + sin2(6'q) sin2(6'p) sin(20p)) 

where Pmax is the largest possible relative proton momentum. This function has been chosen 
to represent a possible, typical dependence of on the final-state kinematics. 

For a large number of arguments x, corresponding to random, phase-space-distributed 
break-up events, the exact value of the function is compared to the value returned by the 
interpolation routine, acting on the grid values. The distribution of differences is shown 
in Fig. ini The spread of values around zero arises because the interpolation assumes that 
the function is linear within a grid cell. This effect causes a systematic shift between the 
interpolated and the true values. For this reason, the FWHM of 0.003 observed in Fig. 1^1 
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should be compared to the value of the function. When averaging over regions of phase 
space this effect washes out further, becoming clearly insignificant indicating that the mesh 
size of the grid used here is sufficiently small. 
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